/*
 * polymer.h
 *
 *  Created on: Oct 3, 2012
 *      Author: hammer
 */

class Polymer {
  int N;
  int DIM;
  vector< vector<float> > positions;
  vector<float> add(vector<float>,vector<float>);
  int increase_i(int i,int j,int l,float sphere_size);
  int decrease_j(int i,int j,int l,float sphere_size);
public:
  Polymer();
  Polymer(int n,int d);
  void add_monomer(vector<float>);
  void print();
  void print(ofstream& out_file);
  int get_length();
  float R_EE();
  float distance(int i,int j);
  vector<float> get_pos(int n);
  void set_pos(int n,vector<float> pos);
  void set_start(vector<float> pos);
  void split(Polymer& poly1,Polymer& poly2,int n);
  void make_random(float step_size);
  Polymer operator+ (Polymer poly1);
  bool check_joining(int join_point,float sphere_size);
  void rotate(float ang);
  void rotate(float ang1,float ang2,float ang3);
  bool pivot(float sphere_size);
};

vector<float> Polymer::add(vector<float> vec1,vector<float> vec2) {
  int i;
  vector<float> vec;
  vec.resize(DIM);
  for (i=0;i<DIM;i++) {
    vec[i]=vec1[i]+vec2[i];
  }
  return vec;
}

Polymer::Polymer() {
  N=0;
  DIM=0;
}

Polymer::Polymer(int n,int d) {
  int i;
  N=n;
  DIM=d;
  positions.resize(N);
  for (i=0;i<N;i++) {
    positions[i].resize(DIM,0);
  }
}


void Polymer::add_monomer(vector<float> position) {
  if (DIM==0) DIM=position.size();
  positions.push_back(position);
  N++;
}

void Polymer::print() {
  int i,j;
  for (i=0;i<N;i++) {
    cout << i <<"\t";
    for (j=0;j<DIM;j++) {
      cout << positions[i][j] <<"\t";
    }
    cout<<"\n";
  }
  cout << "\n";
}

int Polymer::get_length() {
  return N;
}

vector<float> Polymer::get_pos(int n) {
  return positions[n];
}

void Polymer::set_pos(int n,vector<float> pos) {
  positions[n]=pos;
}

void Polymer::set_start(vector<float> pos) {
  if (N==0) add_monomer(pos);
  else {
    int i;
    vector<float> pos1;
    vector<float> diff;
    diff.resize(DIM);
    pos1=get_pos(0);
    for (i=0;i<DIM;i++) {
    	diff[i]=pos[i]-pos1[i];
    }
    for (i=0;i<N;i++) {
      set_pos(i,add(diff,get_pos(i)));
    }
  }
}

void Polymer::make_random(float step_size) {
  int i,j;
  float r;
  bool check;
  vector<float> step;
  step.resize(DIM);
  for (i=1;i<N;i++) {
    check=0;
    while (check==0) {
    	r=0;
    	for (j=0;j<DIM;j++) {
    		step[j]=ran2(&idum)*2-1;
    		r=r+step[j]*step[j];
    	}
    	r=sqrt(r);
    	if (r<1) check=1;
    }
    for (j=0;j<DIM;j++) {
      step[j]=step[j]/r*step_size;
    }
    set_pos(i,add(step,get_pos(i-1)));
  }
}

Polymer Polymer::operator+ (Polymer poly1) {
	int i,n1,n2;
	n1=get_length();
	n2=poly1.get_length();
	poly1.set_start(get_pos(n1-1));
	Polymer temp_poly;
	for (i=0;i<n1;i++) {
		temp_poly.add_monomer(get_pos(i));
	}
	for (i=0;i<n2-1;i++) {
		temp_poly.add_monomer(poly1.get_pos(i+1));
	}
	return temp_poly;
}

void Polymer::print(ofstream& out_file) {
	int i,j;
	for (i=0;i<N;i++) {
		out_file << positions[i][0];
		for (j=1;j<DIM;j++) {
			out_file << "\t" << positions[i][j];
		}
		out_file << "\n";
	}
}

float Polymer::R_EE() {
	int i;
	float r=0;
	for (i=0;i<DIM;i++) {
		r=r+(positions[0][i]-positions[N-1][i])*(positions[0][i]-positions[N-1][i]);
	}
	return r;
}

float Polymer::distance(int i,int j) {
	int k;
	float r=0;
	vector<float> pos1,pos2;
	pos1=get_pos(i);
	pos2=get_pos(j);
	for (k=0;k<DIM;k++) {
		r=r+(pos1[k]-pos2[k])*(pos1[k]-pos2[k]);
	}
	return sqrt(r);
}

bool Polymer::check_joining(int join_point,float sphere_size) {
	float Ri,Rj;
	int i,j,n,count;
	count=0;
	i=join_point+1;
	j=join_point-1;
	while ((j>-1) || (i<N)) {
		count++;
		//cout << "move " << count << ":   i: " << i << ", j: " << j << "\n";
		if (i<N) Ri=distance(i,join_point);
		else Ri=distance(N-1,join_point);
		if (j>-1) Rj=distance(j,join_point);
		else Rj=distance(0,join_point);
		//cout<< "Ri: " << Ri << " Rj: " << Rj << "\n";
		if (((Ri<=Rj) && (i<N)) || ((j<0) && (i<N))) {
			//cout << "increasing i\n";
			n=increase_i(i,j,join_point,sphere_size);
			if (n==-123456789) return 0;
			else i=n;
		}
		else if (((Ri>Rj) && (j>-1)) || ((i>N-1) && (j>-1))) {
			//cout << "decreasing j\n";
			n=decrease_j(i,j,join_point,sphere_size);
			if (n==-123456789) return 0;
			else j=n;
		}
		else cout << "problem";
	}
	//cout << "i: " << i << ", j: " << j << "\n";
	return 1;
}

int Polymer::increase_i(int i,int j,int l,float sphere_size) {
	int j1,int_d,s;
	int b=N-1;
	float step_size=distance(0,1);
	float d;
	j1=l-1;
	if (j<0) j=0;
	while (j1>j-1) {
		d=distance(j1,i);
		if (d<sphere_size*2) {
			//cout << "contact!\nj1= " << j1 << ", i= " << i << ", r= " << d << "\n";
			return -123456789;
		}
		d=d/step_size;
		int_d=floor(d);
		if (int_d>b) s=int_d-b;
		else s=int_d/2;
		if (int_d-s<b) b=int_d-s;
		if (b==0) b=1;
		j1=j1-1-s;
	}
	//cout <<"b= "<< b << ", int_d= " << int_d << ", s= " << s<< "\n";
	return i+b;
}

int Polymer::decrease_j(int i,int j,int l,float sphere_size) {
	int i1,int_d,s;
	int b=N-1;
	float step_size=distance(0,1);
	float d;
	i1=l+1;
	if (i>N-1) i=N-1;
	//cout << "i1: ";
	while (i1<i+1) {
		d=distance(i1,j);
		//cout <<i1 << " ";
		if (d<sphere_size*2) {
			//cout << "contact!\ni1= " << i1 << ", j= " << j << ", r= " << d << "\n";
			return -123456789;
		}
		d=d/step_size;
		int_d=floor(d);
		if (int_d>b) s=int_d-b;
		else s=int_d/2;
		if (int_d-s<b) b=int_d-s;
		if (b==0) b=1;
		i1=i1+1+s;
	}
	//cout <<"\n";
	return j-b;
}

void Polymer::split(Polymer& poly1,Polymer&poly2,int n) {
	int i;
	for (i=0;i<n+1;i++) {
		poly1.add_monomer(get_pos(i));
	}
	for (i=n;i<N;i++) {
		poly2.add_monomer(get_pos(i));
	}
}

void Polymer::rotate(float ang) {
	Fmatrix mat(DIM,DIM,0);
	int i;
	mat.set(0,0,cos(ang));
	mat.set(0,1,-sin(ang));
	mat.set(1,0,sin(ang));
	mat.set(1,1,cos(ang));
	for (i=0;i<N;i++) {
		positions[i]=mat*positions[i];
	}
}

void Polymer::rotate(float ang1,float ang2,float ang3) {
	int i;
	Fmatrix Ax(DIM,DIM,0);
	Fmatrix Ay(DIM,DIM,0);
	Fmatrix Az(DIM,DIM,0);
	Fmatrix A(DIM,DIM,0);
	Ax.set(0,0,1);
	Ax.set(1,1,cos(ang1));
	Ax.set(2,2,cos(ang1));
	Ax.set(1,2,sin(ang1));
	Ax.set(2,1,-sin(ang1));
	Ay.set(1,1,1);
	Ay.set(0,0,cos(ang2));
	Ay.set(2,2,cos(ang2));
	Ay.set(0,2,-sin(ang2));
	Ay.set(2,0,sin(ang2));
	Az.set(2,2,1);
	Az.set(0,0,cos(ang3));
	Az.set(1,1,cos(ang3));
	Az.set(0,1,sin(ang3));
	Az.set(1,0,-sin(ang3));
	A=Az*Ay*Ax;
	for (i=0;i<N;i++) {
		positions[i]=A*positions[i];
	}
}

bool Polymer::pivot(float sphere_size) {
	int i;
	int n=ran2_int(&idum,0,N-1);
	Polymer poly1;
	Polymer poly2;
	Polymer poly3;
	split(poly1,poly2,n);
	if (DIM==2) {
		float ang=ran2(&idum)*2*PI;
		poly2.rotate(ang);
		poly3=poly1+poly2;
		if (poly3.check_joining(n,sphere_size)){
			for (i=n;i<N;i++) {
				set_pos(i,poly3.get_pos(i));
			}
			return 1;
		}
		return 0;
	}
	else if (DIM==3) {
		float ang1=ran2(&idum)*2*PI;
		float ang2=ran2(&idum)*2*PI;
		float ang3=ran2(&idum)*2*PI;
		poly2.rotate(ang1,ang2,ang3);
		poly3=poly1+poly2;
		if (poly3.check_joining(n,sphere_size)) {
			for (i=n;i<N;i++) {
				set_pos(i,poly3.get_pos(i));
			}
			return 1;
		}
		else return 0;
	}
	return 0;
}
